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The stochastic Gross-Pitaevskii equation represents a versatile approach for studying the dynamics 
of trapped degenerate ultracold Bose gases in the presence of large phase and density fluctuations. 
Following a brief review of the original formulation of Stoof, which highlights the benefits of this 
approach and its relation to alternative theories, we present selected applications for the dynamics 
... of effectively one-dimensional systems, and briefly discuss the generalization to two-dimensional 

On ' systems, highlighting certain potential pitfalls in their numerical implementations. 
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I. INTRODUCTION 



The ultracold atomic samples now routinely produced in many laboratories worldwide, provide excellent means 
by which to investigate many-body quantum systems on a macroscopic scale. From a theoretical perspective, these 
Jj^ , systems are appealing as they support an array of rich dynamics, in many cases readily amenable to experimental 
I ' observation. With this close relation between theory and experiment P, 0] , many such experimental features have 
been subject to analytical and numerical investigations: much attention has been paid to numerical analyses based 
upon the Gross-Pitaevskii equation (GPE) in particular 0], which has been shown to successfully describe condensate 
dynamics in the limit when the many-body wavefunction represents the condensate alone, thus corresponding to the 
low-temperature regime. 

^ ' Most experiments on trapped ultracold gases feature, to a greater or lesser extent, a depletion of the condensate due 
to the effect of finite temperature. While in three dimensional systems, finite temperature effects play an increasingly 
^ important role as the critical temperature is approached, for low dimensional systems, thermal effects harbour greater 
O ! importance still, due to the existence of enhanced phase fluctuations within systems under tight confinement in one 
, ^/ or two directions Moreover, the dynamics of the thermal component can play a significant role, for example 

in the case of an ultracold gas under strong external perturbations J, 8, 9, 10]. 
^sj I A number of finite-temperature techniques have been developed to incorporate the effects of temperature into a 
J> ■ description of trapped Bose gases, as recently summarised in [HI, [13, IH, 13 , which feature a detailed comparison of 
the relative weaknesses and merits of such approaches. In brief, various theories have been put forward and applied 
to the experiments, as discussed below: In mean-field type approaches, the condensate is treated as a classical object, 
separated from the effect of the thermal cloud: various such perturbative approaches have been formulated to second 
order in the effective interaction strength [3, [H, [13, [III , with the common feature of obtaining a set of coupled 
dynamical equations for the condensate and the non-condensate, which interact both via mean fields and by particle 
' exchange. The condensate obeys an appropriately generalised GPE, which, in the most commonly adapted version 
I following from early work of Kirkpatrick and Dorfman [l^, [50, [2l| , is coupled to a quantum Boltzmann equation for the 
• • , thermal cloud; this model has been extensively used by Zaremba, Nikuni, Griffin and co-workers, who introduced the 
_ ^ ' terminology 'ZNG' [15] . Although such approaches work remarkably well in predicting a large range of experimentally- 
, relevant phenomena, they do not fully describe regimes where fluctuations in the condensate phase are large, as is 
^ • appropriate for example in low-dimensional systems [1,[^,[1]. To incorporate such effects, generalised zero temperature 
\ [23] and finite temperature [2^ modified methods were developed, to describe the static properties of such systems 
in an ab initio manner. 

However, in order to go beyond mean field effects in a dynamical manner, various alternative approaches have 
been developed and implemented. In first instance, the time-dependent mean field approach mentioned above was 
reformulated in a somewhat more justified number-conserving manner [2^ . following from early work ^2^, [2^ . [30| : 
although this has led to good agreement with experiments in at least one case [31] , a full dynamical model of such 
theories is still lacking. 

Classical field approaches are based upon the observation that [H, [s^], as an equation for a classical field, the GPE 
should be able to describe the classical aspects of an ultracold trapped Bose gas. This includes all highly occupied 
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modes, for which it is possible to approximate the field operators for the creation/annihilation of particles, by c- 
number amplitudes [s^l- Consistency requires only the highly occupied modes to be propagated this way, which can 
be achieved by introducing a projector to the GPE [sH]- Temperature effects can be incorporated through the initial 
conditions, which are typically set to some non-equilibrium distribution, and subsequently allowed to equilibrate at 
some temperature 34] , via the (projected-)GPE, under the restriction of fixed particle number and energy. 

A variant of this technique is known as the Truncated Wigner method, originally developed in the context of 
quantum optics and first used for Bose gases in ; a key element of this method is the approximate incorporation of 
quantum effects, through the addition of a prescribed amount of quantum noise to the wavefunction initial conditions. 
The numerical procedure which results, is to propagate a set of stochastic matter fields, which serve to sample the 
Wigner distribution under this approximation, via an 'appropriate' equation of motion for the Wigner distribution. 
The terminology 'truncated' Wigner arises from the fact that the propagation is carried out in an approximate manner, 
via the (projected) GPE. A limitation of this approach is its inability to accurately describe thermal effects, although 
certain remedies of limited applicability have been proposed [s^l- In the latter two approaches, a temperature can 
only be assigned at the end of the simulations, with the thermalisation arising respectively due to ergodicity [H], 
while simultaneously exhibiting an unphysical heating. 

To gain better control over the temperature to which the system equilibrates, requires detailed consideration of 
the coupling between such modes, and the higher lying, thermal region in our system. Within such schemes, an 
appropriate generalisation has been formulated in f 38|. Earlier, rather distinct, yet analogous in spirit, means of 
achieving this was provided independently by Stoof |39l. [40l . l4ll . |4^ . [i^ . and by Gardiner, Davis and co-authors 
HI, EE S, S^j [13| , using very different theoretical approaches. 

This led to the generalisation of the GPE to a Langevin equation [4l|, often referred to as the Stochastic-Gross- 
Pitaevskii equation (SGPE) [l^, a terminology that we also adopt here; these two independent formulations have 
many formal similarities and are closely related, up to some subtle difference discussed in [l^, EE Ezl- The SGPE 
describes the stochastic evolution of the condensate plus highly occupied, low-lying thermal modes, under the influence 
of coherent and incoherent interactions with the particles occupying higher energy modes. 

Given the diverse range of approaches currently available, our aim in this manuscript is to highlight the key 
elements of this latter approach, paying particular attention on some of its appealing features via somewhat idealised 
experimental scenarios. We start our discussion with a brief overview of one formulation of this approach, which has 
already attracted a relatively large literature [U, E3, El, 113, [HI, [H, [5^ , with related work discussed in |47l . [s^ . [ssl . 



II. METHODOLOGY 



A. The Stochastic Gross-Pitaevskii equation 



The stochastic Gross-Pitaevskii equation is a non-equilibrium microscopic theory to describe the evolution of an 
ensemble of ultracold atoms in contact with a thermal cloud. Using the Keldysh non-equilibrium formalism, Stoof 
derived a Fokker-Planck equation for the the system in (i^l ; such an equation is equivalent to the following Langevin 
equation 41|, 
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where here <l>(r; t) is the order parameter for the lowest, coherent system modes and the higher modes are assumed 
to obey their own coupled dynamics, being in general governed by a quantum Boltzmann equation as in (l^. [isl [l6|. 
In addition to the usual Gross-Pitaevskii terms for the kinetic energy, trapping potential (Vcxt(r)) and nonlinear 
interaction term (where g = ATrh^a^d/'m, parameterised by the s-wave scattering length 03^), this equation features 
two additional contributions: firstly, -R(r; t) denotes a dissipative term which represents the coupling between the 
condensate and a thermal reservoir at a temperature T, and accounts for the transfer of particles between the high 
and low energy modes of the system; the contribution 7/(r; t) is a noise term which accounts for the random nature 
of incoherent collisions within the system and is somewhat analogous to the random particle jitter in Brownian-type 
motion (see, e.g. [l^). Finally, /x denotes an effective chemical potential. The above noise contribution is characterized 
by Gaussian correlations of the form 
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where (...) denotes averaging over different realisations of the noise (and it is generally understood that the results 
of simulations will undergo appropriate averaging, although single-runs can also offer qualitative results as discussed 
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later). The amplitude of the noise is in general position and time dependent, and is governed by the so-called Keldysh 
self-energy S^(r, t). For a thermal cloud which is sufficiently close to equilibrium, this takes the form [4l| 
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Although the occupation numbers, Ni, and therefore also energies Si are in general time-dependent (their values 
governed by a quantum Boltzmann equation and self-consistent Hartree-Fock energies [l^ EH), current numerical 
implementations only focus on near-equilibrium situations, for which the thermal cloud is treated as static; in this case, 
Ni can be approximated by the usual Bose-Einstein distributions Ni = [exp(/3(ei — y))) — 1]~ , where (i — l/ksT 
and Pi labels the momentum of a particle in the ith single particle energy level. 

As one might expect, the strength of the dissipative effects arising from dynamical particle exchange between the 
two subsystems is 'balanced' by the magnitude of fluctuations present, a relation encapsulated by the fluctuation- 
dissipation relation for the system, upon noting explicitly the dependence of these quantities on the system energy 
£c, namely iR{r;ec) — — (l/2)?il]^(r; Cc) [1 -I- 2N{£c)]~^ . In order to make the solution of ([1]) numerically tractable, 
in first instance one can allow the system to equilibrate to a classical distribution by noting that for large oc- 
cupation numbers, [1 + 2N(ec)]~^ ~ l/2/3(ec — A*)- Since the condensate energy is still an operator of the form 
Ec — {—h^y^ /2m -f 14xt(r) + 51^(1"; t)P), the Langevin equation takes the somewhat simplified form 
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The very first numerical implementation of this SGPE in the context of ultracold Bose gases was undertaken in 
2001 by Stoof and Bijlsma [41(| to model the reversible growth of a condensate, as seen in the MIT experiment of 
Stamper-Kurn et al. [58|]: in this experiment, a dimple potential was introduced in a reversible manner to the centre of 
a harmonic trap, thereby inducing localised phase space compression, and leading to a (periodic) crossing through the 
phase transition. Good qualitative agreement was demonstrated between the numerical results (4lj | and the observed 
lagging behind of the condensate growth, relative to the addition of the dimple trap, found in the experiment. The 
method, which is also amenable to analytic variational calculations of stochastic dynamics [13, HH, has since been 
applied (in collaboration with one of us, NPP) to numerical studies of coherence [i^, H^, atom laser [sO], and atom 
chip dynamics '53|. It should be noted that related work has been recently undertaken by Gardiner, Davis and 
co-workers i46|, [13, [H HE HI ■ 

The aim of this paper is to briefly review the main concepts and illustrate the versatility of the SGPE formalism for 
describing systems in which fluctuations are of great importance; for numerical convenience, most of our discussion 
is restricted to effectively one-dimensional systems, for which we also consider the reduction of the SGPE to simpler 
theories, which incorporate finite temperature effects, to differing levels of approximation. The extension of this ap- 
proach to two-dimensional systems is also briefly reviewed. Although, by construction, the results of such simulations 
are best interpreted after numerical averaging over the different noisy initial conditions, we shall see that important 
information is also contained in single runs, as also noted in [sgI. [STj. 



III. ONE-DIMENSIONAL APPLICATIONS 



A. Spatiotemporal Evolution of Coherence in a One-dimensional Bose Gas 

1. Non-equilibrium Coherence 

Low dimensional systems exhibit a richer phase diagram than their three-dimensional counterparts: in the case of 
a purely Id homogeneous system, long wavelength fluctuations in the phase [5^, prevent the onset of Bose-Einstein 
condensation, even down to zero temperature. For a trapped, one-dimensional gas, however, the long-wavelength 
physics is altered due to the low energy cut-off imposed by the harmonic trapping potential, and the quantum 
degenerate regime is instead described by two characteristic temperatures, and T^, below which density and phase 
fluctuations, respectively, become suppressed (25j . The intermediate regime [6^ . is found to be one in which density 
fluctuations no longer domina te, y et fluctuations in the phase remain; this is the regime in which the system is said 
to contain a quasi-condensate [5^ . 
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Figure 1: Snapshots of the density and first-order correlation function during the growth of a one- dimensional Na quasi- 
condensate. The simulations were performed for an effective chemical potential fj. = 30-Eosc, where i?oac = huJz is the harmonic 
oscillator energy, with ujz = 2n x 30Hz at a temperature T = lOOnK for a final quasi-condensate atom number A'^ ~ 20000. 
The effective ID interaction strength gid is obtained by averaging over transverse Gaussian profiles of width li_ — ^/ft/mojx, 
where lu± — 2n x 120H[z. For the chosen somewhat artificial parameters, the phase degeneracy temperature is relatively high, 
slightly in excess of lOOnK. Displayed snapshots were taken at t — 0.05teq, t = 0.15te<j, t = O.Steq and t = OAbteq, where teq 
denotes an approximate time for the system to reach dynamical equilibrium. 



Under tight transverse confinement, the description of the dynamics of highly elongated three-dimensional systems 
effectively reduces to consideration of a one-dimensional equation, as dynamics in the transverse direction becomes 
restricted to the lowest mode. As this is the case in many current experiments with atom chips it is crucial to 
fully understand and characterize such fluctuations. The SGPE is well suited to investigations of phase-fluctuating, 
or quasi-condensates, as fluctuations about the mean fleld are implicitly retained within the order parameter, thus 
giving access to information on coherence, available through the evaluation of various temporal and spatial numerical 
auto-correlation functions. 

A study of the lowest order correlation functions of the gas reveals crucial information about the system. In 
particular, information about the phase of the system is contained in the flrst order correlation function, with the 
next higher order correlation containing crucial details about density fluctuations. Such quantities can be directly 
measured in interferometric 'juggling' experiments (2^. [23l. [2^ . and can also be easily obtained numerically within this 
formalism. More speciflcally, the renormalised first order correlation function at time t can be obtained in 'symmetric 
form' about a point zq, via the relation 

(1). _^ ,^ {<^*{zo-z/2,t)^zo + z/2,t)) 

' (-0 - ./2, zo + z/2; t) . ^|^(,„_,/2,,)|.)(|<,(,„ + ,/2,,)|2) • (4) 

We start our discussion by the basic demonstration of the stochastic method, investigating the nonequilibrium 
evolution of the density and first-order correlation function during the growth of a one-dimensional quasi-condcnsate 
from a static thermal cloud. Fig. [T] shows snapshots of the time evolution of both the density (top row) and coherence 
(bottom row) towards their equilibrium values, for a gas in contact with a heat bath at a fixed temperature. As 
the atomic system is cooled and heads towards the new equilibrium, the spatial extent of the correlation function 
increases as phase coherence is established, with the final shape attained dictated by the temperature, for a given 
particle species, number and trap geometry. In this and all subsequent figures on coherence, distances are scaled to the 
effective quasi-condensate spatial extent Rtf{T) at equilibrium. This quantity takes account of the varying system 
size at different temperatures due to quasi-condensate depletion, and reduces at T = to the 'usual' Thomas-Fermi 
solution Rtf{T = 0)^ muji. There are two ways to obtain Rtf{T), as discussed in [s^: firstly, by considering 

the modified low-dimensional theory of Andersen et al. or equivalently by using the curvature of the density 

profiles to obtain a corresponding quasi-condensate size. The particular example considered here corresponds to the 
case where the coherence length of the system is comparable to the system size, with T w 0.9T^. The observed growth 
in the densities is a direct consequence of the transfer of particles from the heat bath to the system, and the noise is 
crucial in order to seed the growth process [iO, |4lj . 
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Figure 2: Coherence growth at the trap centre for the temperatures T = Q.2T^ and T = 0.9r^. The upper row shows the 
coherence evaluated at a time early in the growth, t = O.lSteij, for each temperature, whereas the plots in the lower row display 
the coherence at the equilibrium time, t — teq. Solid (noisy) curves correspond to numerical results obtained after suitable 
averaging over a few thousand runs, whereas dashed lines indicate the 'optimum' fit using the function f(z) given in the text. 
Other parameters are as in Fig. 1. 



For temperatures such that T > T^, the first-order correlation function has been found experimen tally to exhibit 
exponential decay, whereas for T < the corresponding spatial dependence becomes Gaussian [l^, [IJ. Thus, one 
possible way to quantify the deg ree of coherence and extract a uniformly- varying coherence length is based on a simple 
fitting function of the form [53|, f(z) — exp{—[{z/Lcoh) + C {^1 Lcoh) ]}, in which ^ is a parameter characterizing the 
crossover from exponential to Gaussian behaviour. 

A comparison between the spatial coherence at the trap centre (i.e. for zq — 0) for early and late times, together 
with the results of the fitting procedure, is shown in Fig. [2] for two temperatures, corresponding roughly to the case 
of a 'pure' condensate (left images) and a gas at the crossover from 'pure' to 'quasi'-condensation (right images). 

As is evident from this figure (right images), for the higher temperature results, corresponding to the parameters 
of Fig. [U the equilibrium correlation function decays more rapidly than for the lower temperature limit. Although 
the spatial extent of the system differs between these two cases {Rtf{^-'^T^) = 7.4Z^ vs. i?TF(O.9r0) — 6.6/z), the 
appropriately scaled correlation function plotted in the top row of Fig. [2] exhibit very similar behaviour, unlike the case 
at their respective equilibria (bottom row). Nonetheless, the chosen fitting function seems to appropriately capture 
the main features; importantly, this enables us to study in a universal manner the growth of coherence as the system 
is equilibrating, being driven by the heat bath. 

Numerous experiments [6l|, [gI, [H, [131 have been performed to study the relation between condensate growth and 
evolution of coherence, and a detailed comparison to such experiments is pending (but see also [H,!!!] and references 
therein). For our present illustrative purposes, we simply highlight that the growth of coherence appears to exhibit 
similar behaviour to typical quasi-condensate number growth curves, although these generally follow such curves with 
an additional small lagging time; most notably, once the norm has reached its steady-state value, the coherence is 
still growing, albeit at a very slow rate. 

With this means of quantifying the degree of coherence across a sample, it is then interesting to assess how Lcoh 
varies with time during the cooling process. The analysis of these results for the two different temperatures in the 
context of this somewhat idealised scenario is plotted in Fig. [3l As these simulations are based on a fixed value of /z, 
the final atom numbers are slightly increased in the high temperature case (top right image) ; note also the significantly 
increased relaxation time and observation of the initial spontaneous initiation process prior to the domination of the 
bosonic stimulation effect (top left). 
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Figure 3: Quasi-condensate atom number growth curves for a system equilibrating in contact with a prescribed heat bath (top 
row), versus corresponding temporal evolution of the coherence length Lcoh, obtained to approximately 10% accuracy from 
the fitting function defined in the text (bottom row). Note the differing y-axis scales between the coherence growth for two 
temperatures, with results scaled to the temperature dependent Thomas- Fermi radius Rtf{T). 




z/Rtf(T) z/R^p(T) 

Figure 4: Evolution of the coherence (from left to right within each image) during the growth to equilibrium with the heat bath 
for a temperature T = O.QT^. The correlation function is evaluated at zo — 0.7Rtf{T) (lower left plot) and zq = 1.4:Rtf{T) 
(lower right plot), with zo indicated by the circles on the equilibrium density profile for the chosen temperature (top plot). 
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2. Dependence of Equilibrium Coherence on Position 

We also investigate the position dependence of the correlation function at equilibrium, by comparing the evolution 
of g^^^ at different positions zq, for a fixed temperature corresponding to the high temperature case presented in the 
above results. The results of this investigation are shown in Fig. [Hfor z = O.TRtf and z = IARtf for various times 
during the growth to equilibrium. The final profile attained at this temperature, together with the Thomas-Fermi 
profile for the same number of atoms, is also shown, with the positions considered indicated. The coherence within 
the Thomas-Fermi radius decreases uniformly but at a relatively slow rate as soon as one probes off-centred regions, 
with the effective coherence length decreasing rapidly as soon as one moves beyond the Thomas-Fermi radius (right 
images). 



3. Potential 'Identification' of a Quasi- condensate 



An appealing feature of the stochastic approach that is often not fully appreciated, is that it actually generates 
total density profiles without needing to resort to a special mean-field treatment of the condensate mode. This means 
that, as in the experiments, one does not have an automatic 'visualisation' of the condensate in the system. So, in 
order to analyse these results, one may consider applying bimodal fits precisely as done in the experiments. However, 
an alternative definition has been proposed [13], in which the quasi-condensate density is identified as riQci.z) = 
\/ {\^{z)\^)'^ — {\^{z)\^). This approach has been implemented in [5^, with the density profile contribution for the 
thermal cloud over low-lying modes determined self-consistently via the relation nT(-2^) — (I'^'C-^)!^) — "-qcC-^), within 
the region z < Rtf{T)^ and given simply by ut^z) = n{z) outside of the Thomas-Fermi radius. The identification of 
these quantities enables spatial plots of the temperature dependent quasi-condensate density profiles to be obtained 
in an ab initio manner, as shown for different temperatures in Fig. [5] Such a definition was used more recently to 
distinguish between a condensate and a quasi-condensate [H, [sH] 

Having discussed some basic elements of finite temperature trapped, quasi-one-dimensional Bose gases via the 
SGPE, we now discuss certain limiting cases of the above theory which are commonly used in the literature, and 
which shed more light on the benefits and versatility of the SGPE. 




Figure 5: Quasi-condensate (solid) and thermal (dashed) density profiles obtained from the SGPE by manipulation of density- 
density correlations. Plotted images correspond to the following temperatures T/T^ — 0.1, 0.9, 1.8, and 2.9 (top left to bottom 
right). 
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B. Reduction to Alternative Theories 



We have already argued that the two important features of the SGPE are that it includes the correct damping due 
to its dynamical coupling to the heat bath, and that it extends beyond mean field effects, by including a stochastic 
noise term. Removing either (or both) of these elements leads to simpler theories [12| , which can nonetheless be useful 
in certain regimes, as discussed below. 



1. Removing Contact with the Heat Bath 

In first instance, we can consider what happens when the coupling to the heat bath is abruptly (and perhaps 
somewhat artificially) removed. We recall that upon beginning a simulation with the SGPE, the effect of the iR{r;t) 
term in ([T]) is to pump particles from the heat bath, and populate the low lying modes to their equilibrium occupations, 
in accordance with the classical equilibrium for a given set of trap frequencies, chemical potential and temperature. 
On reaching a dynamical equilibrium, the scattering rate oc iR{r; t) becomes on average equal to zero, and there is no 
net particle exchange with the heat bath, although this continuous exchange can cause additional damping. However, 
any externally imposed perturbations of the system about this set of system parameters, would typically lead to the 
introduction or loss of particles, upon further evolution with ([1]). 

In general, the SGPE is coupled to a quantum Boltzmann equation for the population of the higher-lying modes, 
and these two equations should be solved self-consistently. If this is not done, then modifying for example the trap 
geometry will lead to a change in the norm of the SGPE. If we are only considering the evolution of the low-lying 
modes, in order to keep the atom number within our system fixed, one could consider 'fixing' the norm in the numerics 
by removing the coupling to the heat bath. A better alternative might perhaps be to change the effective chemical 
potential in the SGPE to the final value corresponding to the same atom number in the new trap, although the 
observed dynamics will typically depend on how quickly the chemical potential is adjusted to its final value. 

One procedure is thus to propagate initially the full stochastic solution in order to obtain the desired equilibrium 
state, before switching off the noise and damping terms of ([1]) , and subsequently propagating the noisy initial condition 
via the usual GPE which enforces particle number conservation. This corresponds physically to the removal of contact 
with the heat bath upon reaching equilibrium. 

As a method, this procedure has strong parallels with the Truncated Wigner method, and a comparative study 
of the two is under way curre ntly [6 81. The main difference arises in the fact that in the usual truncated Wigner 
implementations (see, e.g. [36ll69ll70l| or [T^ for a more extensive list), only quantum (as opposed to thermal) noise is 
typically added in the initial conditions (but see also [l^jlzii), whereas the SGPE in general contains both quantum 
and thermal fluctuations [1^ . Although existing numerical implementations of the SGPE focus on the thermal effects 
(due to the 'classical' approximation mentioned earlier), the important advantage of this approach is that there is 
an explicit fluctuation-dissipation relation which guarantees convergence to the correct classical equilibrium at any 
pre-determined temperature; in other words, the initial thermal state for the unperturbed system can be correctly 
assigned. The validity of conventional truncated Wigner implementations is on the other hand limited to very low 
temperatures [s^l and occupation numbers for which quantum fluctuations have a dominant role - see [1^ for a more 
detailed discussion. 

1 n ' ' ' ' ' ' ' 1 
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Figure 6: Suppression of Shockwaves induced upon the introduction of a deep dimple trap to a one-dimensional ultracold 
Bose gas of around 3500 ^''Na atoms, as manifested in the temporal evolution (in units of the inverse trap frequency) of the 
corresponding damped central density oscillations obtained after numerical averaging, and scaled to the peak density for the 
given system. (Details of other parameters used in these simulations can be found in [s^.) Contact with the heat bath was 
removed prior to the introduction of the dimple trap to maintain number conservation. 
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This finite temperature approach based on the removal of the heat bath once the correct equihbrium has been 
obtained via the SGPE was used to probe the dynamics of an initially equilibrated thermal stationary state following 
the addition of a deep Gaussian dimple trap [s^] ■ A large and sudden perturbation is found to lead to the appearance 
of shock waves in the system, and in the averaged profiles typically studied, these manifest themselves in terms of 
an initial increase in the central density, followed by periodic damped decrease-increase cycles, as shown in Fig|6l 
this can be interpreted as direct evidence of the temperature-induced suppression of the initial Shockwave amplitude 
induced by the addition of the new trap. It should be noted here that the introduction of the correct thermal noise 
in the initial conditions is crucial to the generation of the damping which would not be present in simulations of the 
ordinary GPE starting from a stationary condensed state. 



2. A Microscopically-based Dissipative GPE 

An even simpler limiting case of the SGPE is based on removing only the noise term of (IT|); then, the SGPE reduces 
to a 'dissipative GPE', as noted for example in [13, [Z^- The form of the damping is given by the self energy, and as 
such is spatially dependent, consistent with the spatially inhomogeneous background of the harmonic trap. Such a 
model has been used to study vortex and soliton dynamics leading to a number of interesting predictions. However, 
the origin of such models can be traced to rather crude approximations in the more advanced microscopic models, 
such as the SGPE, with most su ch p ublications justifying this equation phenomenologically, and actually using a 
constant value for the damping 7 [73l. ItI. ItsI. [tgI] . 




Time / »' 



Figure 7: Soliton dynamics within a one-dimensional gas oi N ^ 20000 ^^Na atoms, for the zero temperature (GPE) case, 
plus at temperatures of lOOnK and 200nK. The top row plots show density snapshots at three times during the dynamics, 
t — 24a;J^, t — 37.2uj^^ and, t — 60io^^; shown in each figure are the results for all three temperature cases at that time (zero 
temperature - blue, dot-dashed line; lOOnK - green, dashed line; 200nK - red, solid line). The spacetime plots are ordered 
downwards in order of increasing temperature, for the highest temperature case the soliton can be seen to decay rather rapidly. 

The dissipative equation which results from the SGPE in general takes, within the 'classical approximation', the 
form 



2m 



$(r;t) , 



(5) 
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where the damping term (for a system made up of a fixed total number of atoms) is given by 

7(r,T,i) = *fftS^(r,t), (6) 

ahhough our discussion so far has additionally ignored the time-dependence for numerical simplicity. 

It is very easy to see that the resulting evolution is no longer norm conserving, and it therefore becomes nessessary 
to renormalize the wavefunction as each time step to maintain a constant particle number. 

For illustrative purposes, we apply this approach here to model the oscillations of an 'ideally generated' dark 
soliton within a one-dimensional condensate. The dynamics for three identical initial soliton solutions are shown at 
different temperature in Fig. [7] (top images), with each image showing the position of the dark soliton for three 
different temperatures, corresponding to a different evolution time. To best illustrate the different temperature- 
dependent decay rates, we give corresponding space-time plots of the density along the z-axis with such damping 
already observed in experiments 77]. 

Although consideration of the ID case enables us to bring out the main physical features of the SGPE model, 
comparison to experiments requires us to consider its properties in higher-dimensional systems, and we conclude our 
presentation by a simple illustration in two-dimensional geometries, and some related comments. 



IV. EXTENSION TO TWO DIMENSIONS 



As in one dimension, for a uniform two dimensional Bose gas, thermal fluctuations remain even at very low temper- 
atures, preventing the onset of Bose-Einstein condensation |7 8| . l79l . For interacting trapped two dimensional gases, 
however, quasi-condensation is possible for low temperatures [80|, with a Berezinskii-Kosterlitz-Thouless phase tran- 
sition to a superfluid state associated with the the pairing of vortices of opposite circulation beneath some critical 
temperature. Recently, this has spawned much interest in the physics of two-dimensional trapped ultracold systems, 
particularly in assessing the nature of the phase transition at the critical point [1, [Sll, [H, H^l . 



A. Spontaneous Processes - Vortex Nucleation 



Phase transitions are st rong ly associated with such topological defects, the appearance of which are predicted by the 



Kibble-Zurek mechanism [8J, |85| . An appealing feature of the stochastic approach is that in including fluctuations 



about the mean field, excitations in the form of spontaneous vortices can be seen during the growth process, as 
investigated recently using an alternative formulation of the SGPE [s^- Our discussion so far has been restricted to 
analysing predictions of the SGPE based on numerical averaging over different initial noisy conditions. However, the 
monitoring of the evolution within a single run can also provide important information on the system parameters, 
which may otherwise be lost through averaging - e.g. information related to spontaneous events, such as topological 
defect formation (86j . 

Within the numerical SGPE applied here, condensate formation from the static thermal cloud is initiated for a 
positive value of the chemical potential in A rapid quench can be simulated by an instantaneous change in the 
chemical potential, leading to a rapid condensate growth. Within simulations, the appearance of a (random) number 
of vortices is common, although their details vary from run to run. The appearance of these actually becomes 'washed 
out' following the averaging process, as, due to the random nature of the nucleation of vortices, the positions at 
which they appear varies between different realisations. Such studies should be contrasted to studies based on the 
GPE, in which a vortex is actually 'artificially' imprinted at a predetermined position; when modelling the ensuing 
dynamics, one often introduces damping by resorting to the dissipative GPE ([5]), with the damping term, ^(x,y), 
which is actually proportional to the two-dimensional self-energy at a particular temperature, typically treated as a 
constant. 

In density plots, vortices appear as dark regions, corresponding to positions at which singularities in the phase of 
the system can be seen. An example of the results from a typical simulation based on the SGPE is shown in Fig. 
|51 In the density and phase snapshots shown, two spontaneously generated vortices (top left) can be seen to come 
together (middle) and then annihilate (right), illustrating the spontaneous dynamics captured within simulations. 
The averaged density profiles included show the result of averaging over 100 runs, and the smooth density plots which 
result. 
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Figure 8: Nonequilibrium density and phase snapshots for the growth of a two dimensional harmonically trapped quasi- 
condensate in the x-y plane, for A'' ~ 60000 ^'^Na atoms, dimensionless coupling constant g2d = V^o-sd/lz ~ 0.042 and trap 
frequencies oj^ = ojy = 27r x 200Hz and tJz = 27r x 4kHz. The single run density profiles (top row) show the spontaneous 
appearance of vortices, as apparent also in the single run phase plots (middle row). Time increases to the right in the plots; 
at the centre of the single run images two vortices can be seen to come together and annihilate. The bottom row profiles are 
density profiles obtained following averaging over 100 single runs, with snap-shots taken at the same times as for the single 
runs. The results show a far smoother density profile is obtained upon averaging, though vortices are now no longer visable. 




0.2 - 



I ' ' ' ' ' ' 1 
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Grid spacing / 1^ 

Figure 9: Dependence of the norm on the numerical grid spacing for the two-dimensional SGPE. The fit shows the dependence 
of the particle number on the high momentum cut-off, A, introduced by the lattice for fixed chemical potential, temperature 
and trap frequencies. This is inversely related to the grid spacing Az, against which the results are plotted. The divergence is 
found to be oc log(A) in the two dimensional case, as shown by the line indicating the fit to a function cx log(7r/A2). 
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B. Numerical Issues - Dependence on Momentum Cutoff 



The numerical means of evolving the system discussed so far make use of the macroscopic occupation of the low 
lying modes, assuming these states to be well described by a classical field. In two dimensions, known divergences 
arise, associated with large momenta, or equivalently small grid spacings, as the continuum limit is approached on the 
lattice. This is due to the high momentum cutoff, A, related to the numerical grid spacing, Az, through the relation 
A = tt/Az. 

We present here preliminary results of an investigation into this dependence for the SGPE as a 'caution' to the 
reader. Here we make no adjustments to parameters to compensate for grid dependence and find the expected 
logarithmic divergence in the momentum cut-off, with the system norm (or number of particles in the classical region) 
increasing with decreasing spatial discretization Az, as shown in Fig. [9l 

In principle such divergences can be removed from the problem through the appropriate renormalization of physical 
quantities from their original ('bare') values. For example in homogeneous systems, the lattice effects on temperature 
or density can be taken into account in an a posteriori manner, as described in (STj . An alternative method is 
to introduce counter-terms to the system Hamiltonian, which has been successfully applied to Langevin equations 
previously in, for example, fssl. [sol [ool [9]| . 



The stochastic Gross-Pitaevskii is emerging as a key approach for the equilibrium and dynamical properties of 
degenerate Bose gases, and is particularly beneficial for describing regimes of large phase and density fiuctuations, 
such as near (or at) the phase transition, or for weakly-interacting low-dimensional Bose gases. Single run results 
contain important (at least qualitative) information on spontaneous processes, such as defect formation characteristics, 
whereas averaged profiles are most suitable for determination of correlation functions, quasi-condensate fractions and 
smooth density profiles. Although these stochastic simulations have obvious advantages, their existing numerical 
implementation have not yet explored such theories to their full potential, with current models typically assuming 
a static (near-equilibrium) thermal cloud, and additionally discarding quantum fiuctuations in simulations. Finally, 
the important issue of the cut-off dependence, already highlighted by other researchers should not be overlooked, 
particularly for dimensions higher than one, and various techniques are available for renormalising the observed 
profiles to their actual values, thereby eliminating (at least to leading order) such dependences. We believe that use 
of the stochastic Gross-Pitaevskii equation, which is also coupled to a quantum Boltzmann equation for the self- 
consistent dynamical treatment of the higher-lying thermal modes, should be able to describe essentially all main 
features of weakly-interacting Bose gases currently pursued experimentally, and therefore look forward to the exciting 
times which lie ahead. 



We are grateful to Henk Stoof for early collaboration and discussion on these ideas and to the UK EPSRC for 
funding. 



[1] Dalfovo F, Giorgini S, Pitaevskii L P, and Stringari S 1999 Rev. Mod. Phys. 71 463 
[2] Bloch I, Dalibard J and Zwerger W 2008 Rev. Mod. Phys. 80 885 

[3] Emergent Nonlinear Phenomena m Bose-Einstein Condensates: Theory and Experiment, Ed. Kevrekidis P G, Frantzeskakis 
D J and R. Carretero-Gonzalez (Springer Verlag, 2008) 

[4] Gorlitz A, Vogels J M, Leanhardt A E, Raman C, Gustavson T L, Abo-Shaeer J R, Chikkatur A P, Gupta S, Inouye S, 
Rosenband T and Ketterle W 2001 Phys. Rev. Lett. 87 130420 

[5] Ott H, Fortagh J, Schlotterbeck G, Grossmann A and Zimmermann C 2001 Phys. Rev. Lett. 87 230401 

[6] Hadzibabic Z, Kriiger P, Cheneau M, Battelier B and Dalibard J 2006 Nature (London) 441 1118 

[7] Jackson B and Zaremba E 2001 Phys. Rev. Lett. 87 100404 

[8] Jackson B and Zaremba E 2002 Phys. Rev. Lett. 89 150402 

[9] Jackson B and Zaremba E 2002 Phys. Rev. Lett. 88 180402 

[10] Jackson B and Zaremba E 2002 Laser Physics 12 93 

[11] Proukakis N P Beyond the Gross- Pitaevskii Mean Field in Ref. [3] [ arXiv:0706.35 4lYl ] 
[12] Proukakis N P and Jackson B 2008 J. Phys. B 41 203002 




V. 



CONCLUSIONS 



Acknowledgments 



13 



[13] Blakie P B, Bradley A S, Davis M J, Ballagh R J and Gardiner C W 2008 Advances in Physics 57 363 
[14] Yukalov V 1 2007 Las. Phys. Lett. 4 632 

[15] Zaremba E, Nikuni T and Griffin A 1999 Low J Temp. Phys. 116 277 
[16] Bijlsma M J, Zaremba E and Stoof H T G 2000 Phys. Rev. A 62 063609 
[17] Walser R, Williams J, Cooper J and Holland M 1999 Phys. Rev. A 59 3878 
[18] Proukakis N P 2001 J. Phys. B 34 4737 

[19] Kirkpatrick T R and Dorfman J R 1983 Phys. Rev. A 28 2576 

[20] Kirkpatrick T R and Dorfman J R 1985 J Low Temp. Phys. 58 301 

[21] Kirkpatrick T R and Dorfman J R 1985 J. Low Temp. Phys. ibid. 58 399 

[22] Richard S, Gerbier F, Thywissen J H, Hugbart M, Bouyer P and Aspect A 2003 Phys. Rev. Lett. 91 010405 
[23] Hellweg D, Gacciapuoti L, Kottke M, Schulte T, Sengstock K, Ertmer W, and Arlt J J 2003 Phys. Rev. Lett. 91 010406 
[24] Gacciapuoti L, Hellweg D, Kottke M, Schulte T, Ertmer W, Arh J J, Sengstock K and Santos L 2003 Phys. Rev. A 68 
053612 

[25] Petrov D S, Shlyapnikov G V, and Walraven J T M 2000 Phys. Rev. Lett. 85 3745 

[26] Andersen J O, Al Khawaja U and Stoof H T G 2002 Phys. Rev. Lett. 88, 070407 

[27] Gardiner S A and Morgan S A 2007 Phys. Rev. A 75 043621 

[28] Girardeau M D and Arnowitt R 1959 Phys. Rev. 113 755 

[29] Gardiner G W 1997 Phys. Rev. A 56 1414 

[30] Gastin Y and Dum R 1998 Phys. Rev. A 57 3008 

[31] Morgan S A, Rusch M, Hutchinson D A W and Burnett K 2003 Phys. Rev. Lett. 91 250403 

[32] Kagan Y and Svistunov B V 1994 Zh. Eksp. Teor. Fiz. 105 353 [1994 Sov. Phys. JETP 78 187] 

[33] Brewczyk M, Gajda M and Rzazewski K 2007 J. Phys. B 40 Rl 

[34] Blakie P B and Davis M J 2005 Phys. Rev. A 72 063608 

[35] Davis M J, Morgan S A and Burnett K 2001 Phys. Rev. Lett. 87 160402 

[36] Steel M J, Olsen M K, Plimak L I, Drummond P D, Tan S M, GoUett M J, Walls D F and Graham R 1998 Phys. Rev. A 
58 4824 

[37] Sinatra A, Lobo G and Gastin Y 2002 J. Phys. B 35 3599 

[38] Davis M J, Ballagh R J and Burnett K 2001 J. Phys. B 34 4487 

[39] Stoof H T G 1997 Phys. Rev. Lett. 78 768 

[40] Stoof H T G 1999 Low J Temp. Phys. 114 11 

[41] Stoof H T G and Bijlsma M J 2001 J. Low Temp. Phys. 124 431 

[42] Duine R A and Stoof H T G 2001 Phys. Rev. A 65 013603 

[43] Stoof H T G in Coherent Atomic Matter Waves, Proceedings of the Les Houches Summer School Session 72, 1999, ed. 

Kaiser R et al. (Springer- Verlag, Berlin 2001). 
[44] Gardiner G W and ZoUer P 2000 Phys. Rev. A 61 033601 
[45] Davis M J, Gardiner G W and Ballagh R J 2000 Phys. Rev. A 62 063608 
[46] Gardiner G W Anglin J R, and Fudge T I A 2002 J. Phys. B 35 1555 
[47] Gardiner G W and Davis M J 2003 J. Phys B 36 4731 

[48] Bisset R N, Davis M J, Simula T P, and Blakie P B 2009 Phys. Rev. A 79 033626 

[49] Al Khawaja U, Andersen J O, N. Proukakis P, and H.T.G. Stoof 2002 Phys. Rev. A 66, 013615 ; ibid. 66, 059902(E) (2002) 
[50] Proukakis N P 2003 Las. Phys. 13 527 

[51] Duine R A, Leurs B W A and Stoof H T G 2004 Phys. Rev. A 69 053623 

[52] Proukakis N P, Schmiedmayer J and Stoof H T G 2006 Phys. Rev. A 73 053603 

[53] Proukakis N P 2006 Phys. Rev. A 74 053617 

[54] Bradley A S, Blakie P B and Gar diner G W 2005 J. Phys. B 38 4259 
[55] Bradley A S, and Gardiner G W, |cond-mat/0602162| 

[56] Bradley A S, Gardiner G W and Davis M J 2008 Phys. Rev. A 77 033616 

[57] Weiler G N, Neely T W , Scherer D R, Bradley A S, Davis M J, Anderson B P 2008 Nature 455 948 
[58] Stamper-Kurn D M, Miesner H-J, Ghikkatur A P, Inouye S, Stenger J and Ketterle W 1998 Phys. Rev. Lett. 81 2194 
[59] Popov V N 1965 Sov. Phys. JETP 20 1185; Popov V N Functional Integrals and Collective Excitations (Cambridge 
University Press, 1987) 

[60] Al Khawaja U, Proukakis N P, Andersen J O, Romans M W J and Stoof H T G 2003 Phys. Rev. A 68 043603 
[61] Miesner H J, D.M. Stamper-Kurn, Andrews M R, Durfee D S, Inouye S and Ketterle W 1998 Science 279 1005 
[62] Kohl M, Davis M J, Gardiner G W, Hansch T W and Esslinger T 2002 Phys. Rev. Lett. 88 080402 

[63] Shvarchuck I, Buggle Gh., Petrov D S, Dieckmann K, Zielonkowski M, Kemmann M, Tiecke T G, von Klitzing W, 

Shlyapnikov G V and Walraven J T M 2002 Phys. Rev. Lett. 89 270404 
[64] Ritter S, Ottl A, Donner T, Bourdel T, Kohl M and Esslinger T 2007 Phys. Rev. Lett. 98 090402 
[65] Hugbart M, Retter J A, Varon A F, Bouyer P, Aspect A and Davis M J 2007 Phys. Rev. A 75 011602(R) 
[66] Davis M J and Gardiner G W 2002 J. Phys. B 35 733 
[67] Prokof'ev N and Svistunov B 2002 Phys. Rev. A 66 043608 
[68] Gockburn S P, Negretti A, Henkel G and Proukakis N P 2009 (in preparation) 
[69] Norrie A A, Ballagh R J and Gardiner G W 2005 Phys. Rev. Lett. 94 040401 
[70] Scott R G, Hutchinson D A W and Gardiner G W 2006 Phys. Rev. A 74 053605 
[71] Ruostekoski J and Isella L 2005 Phys. Rev. Lett. 95 110403 



14 



[72] Penckwitt A A, Ballagh R J and Gardiner C W 2002 Phys. Rev. Lett. 89 260402 

[73] Choi S, Morgan S A and Burnett K 1998 Phys. Rev. A 57 4057 

[74] Tsubota M, Kasarnatsu K and Ueda M 2002 Phys. Rev. A 65 023603 

[75] Kasarnatsu K, Tsubota M and Ueda M 2003 Phys. Rev. A 67 033610 

[76] Proukakis N P, Parker N G, Barenghi C F and Adams G S 2004 Phys. Rev. Lett. 93 130408 

[77] Burger S, Bongs K, Dettmer S, Ertmer W, Sengstock K, Sanpera A, Shlyapnikov G V and Lewenstein M 1999 Phys. Rev. 
Lett. 83, 5198 

[78] Mcrrnin N D and Wagner H 1966 Phys. Rev. Lett. 22 1133 

[79] Hohenberg P C 1967 Phys. Rev. 158 383 

[80] Petrov D S, Holzmann M, and Shlyapnikov G V 2000 Phys. Rev. Lett. 84 2551 

[81] Hadzibabic Z, Kriiger P, Cheneau M, Rath S P and Dalibard J 2008 New J. Phys. 10 045006 

[82] Schweikhard V, Tung S, and Gornell E A 2007 Phys. Rev. Lett. 99 030401 

[83] Simula T P and Blakie P B 2006 Phys. Rev. Lett. 96 020404 

[84] Kibble T W B 1976 J. Phys. A 9 1387 

[85] Zurck W H 1985 Nature (London) 317 505 

[86] Anglin J R and Zurek W H 1999 Phys. Rev. Lett. 83 1707 

[87] Prokof'ev N, Ruebenacker O, Svistunov B 2001 Phys. Rev. Lett. 87 270402 

[88] Borrill J and Gleiser M 1996 Nucl. Phys. B 483 416 

[89] Gagne G and Gleiser M 2000 Phys. Rev. E 61 3483 

[90] Wojtas D, 2006 MSc. Thesis, University of Canterbury, Christchurch 

[91] Lythe G, Bettencourt L and Habib S 1999 Physical Review D 60 105039 



